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Preface 


This semi-annual report describes our progress during the period from February 
1990 to September 1990. Several technical reports and papers have been written and these 
are listed at the end of each task. 

There are two tasks described in this report. Each should be read independently. 
That is, figure and reference numbering is consecutive only within the description of the 
task. As can be expected, the progress reports are very brief and the reader should refer to 
the referenced technical reports for detailed coverage. A total of sixteen technical reports 
have already been submitted and two more are currently being prepared. Also, more than 
18 publications to refereed journals and more than 17 conference papers have resulted from 
this sponsored research. 




TABLE OF CONTENTS 


Page # 

PREFACE i 

TASK : A FINITE ELEMENT CONJUGATE GRADIENT FFT 

METHOD FOR SCATTERING 

Abstract *■ 

Objective * 

Progress 3 

Pre- and Post-Processing Algorithms 3 

3D Algorithms for Bodies of Revolution 5 

Conclusions ® 

Transitions ° 

References ' 

Publications Resulting from this Research 10 

TASK: ANALYTICAL SOLUTIONS WITH GENERALIZED 

IMPEDANCE BOUNDARY CONDITIONS (GIBC) 

Abstract 1 3 

Objective 1 3 

Introduction 1 3 

Dual Integral Equation Formulation 15 

Solution of the Diffracted Field 20 

Diffraction by Thin Single Layer Discontinuous Strips 22 

Modal Decomposition of the Symmetric Slab Fields 25 

Exterior Cross Section Functions 28 

Interior Cross Section Functions 29 

Recasting of the Dual Integral Equation Solution for 
a Material Junction 3 1 



Determination of the Constants 34 

Validation of the Solution 37 

Other Applications of the GIBC/GSTC 37 

Appendix 39 

References 41 

Publications Resulting from this Research 44 


^91 - 1447 


Task Title: 


Investigators : 


A Finite Element Conjugate Gradient 
FFT Method for Scattering 


Jeffery D. Collins, John Zapp, Chang-Yu Hsa and John L. Volakis 


Period Covered: 


February 1990 - September 1990 



ABSTRACT 


An extension of the two-dimensional formulation developed last year is presented 
for a three dimensional body of revolution. With the introduction of a Fourier expansion of 
the vector electric and magnetic fields, a coupled two-dimen ional system is generated and 
solved via the finite-element method. As before, an exact boundary condition is employed 
to terminate the mesh and the FFT is used to evaluate the boundary integrals for low O(n) 
memory demand when an iterative solution algorithm is used. Again, by virtue of the Finite 
element method the algorithm is applicable to structures of arbitrary material composition. 

Several improvements to our two-dimensional algorithm are also described. These 
include (1) modifications for terminating the sh at circular boundaries without distorting 
the convolutionality of the boundary ■. tv. egr "Us, (2) the development of our own non- 
proprietery mesh generation routines for two-dimensional applications, (3) the development 
of preprocessors for interfacing SDRC IDEAS@ with the main algorithm, and (4) the 
development of post-processing algorithms based on the public domain package GRAFIC 
to generate 2D and 3D gray level and color field maps. 


OBJECTIVE 

The objective of this task is to develop innovative techniques and related software 
for scattering by three dimensional composite structures. The proposed analysis is a hybrid 
finite element-boundary integral method formulated to have an O(n) memory demand. This 
low storage is achieved by employing the FFT to evaluate all boundary integrals and 
resorting to an iterative solution algorithm. Particular emphasis in this task is the 
generation of software applicable to airborne vehicles and the validation of these by 
comparison with measured and other reference data. Because the approach is new, a step 
by step development procedure has been proposed over a three-year period. During the 
first year the technique was developed and implemented for two-dimensional composite 
structures. Support software for the two-dimenional analysis such as pre- and post- 
processor routines were developed during the second year and a formulation was also 
developed and implemented for three-dimensional bodies of revolution. Finally, during the 
third year, we will develop, implement, and test the method for arbitrary three dimensional 
structures. 
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BACKGROUND 


Interest in three-dimensional (3-D) methods has increased in recent years, however, 
the associated demands in computation time and storage are often prohibitive for electrically 
large 3-D bodies. Vector and concurrent (i.e. hypercube, connection, etc.) computers are 
beginning to alleviate the first of these demands, but a minimization of the storage 
requirements is essential for treating large structures. 

The traditional Conjugate Gradient Fast Fourier Transform (CGFFT) method [ 1 ]- 
[4] is one such frequency domain solution approach which requires O(n) storage for the 
solution of n equations. This method involves the use of FFTs whose dimension equals 
that of the structure under consideration [5]-[7] and, therefore, demands excessive 
computation time when used in an iterative algorithm. Also, the standard CGFFT requires 
uniform rectangular gridding that unnecessarily includes the impenetrable portions of the 
scatterer. With these issues in mind, a new solution approach is propsed for solving 
scattering problems. The proposed method will be referred to as the Finite Element- 
Conjugate Gradient Fast Fourier Transform (FE-CGFFT) method. 

During last year's effort the FE-CGFFT method was developed for two- 
dimensional scatterers where the finite element mesh was terminated at a rectangular box. 
Inside the box boundaries, Helmholtz equation is solved via the finite element method and 
the boundary constraint is obtained by an appropriate integral equation which implicitly 
satisfies the radiation condition. Along the parallel sides of the box, this integral becomes a 
convolution and is, therefore, amenable to evaluation via the FFT. The dimension of the 
required FFT in this hybrid method is one less than the dimensionality of the structure thus, 
making it attractive for 3-D simulations. Also, because it incorporates the finite element 
method, the FE-CGFFT formulation remains valid regardless of the structure's geometry 
and material composition. 

The proposed method described in the University of Michigan Report 025921-6-T 
(see also [8]) is similar to the moment method version developed by Jin [9]. Jin’s method 
was in turn based on work published in the early 70's by McDonald and Wexler [10] who 
introduced an approach to solve unbounded field problems. The proposed method is also 
similar to other methods (a few of which will be mentioned here), neither of which 
provides a storage reduction comparable to the proposed FE-CGFFT method. The 
unimoment method [11] uses finite elements inside a fictitious circular boundary and an 
eigenfunction expansion to represent the field in the external region. The coefficients of the 
expansion are then determined by enforcing field continuity at the finite element (FE) mesh 
boundary. The coupled finite element-boundary element method [12] uses the finite 
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element method within the boundary and the boundary element method to provide the 
additional constraint at the termination of the mesh. Unlike the proposed method, the 
solution in [12] was accomplished by direct matrix inversion (as in [9]), and the outer mesh 
boundary is not rectangular to take advantage of the FFT for the evaluation of the boundary 
integrals. 


PROGRESS 

The proposed FE-CGFFT formulation was implemented last year (see Figs. 1 and 
2) but as can be expected, the rectangular mesh boundary does not always lead to the most 
efficient formulation, particularly when dealing with structures whose outer boundary is 
not rectangular. Because of this, during this year we developed and implemented a 
formulation which permits mesh termination at circular (see Fig. 3) boundaries for the 2D 
case with the corresponding boundary enclosure being a pillbox for the 3D case (see Fig. 

4). As before, these boundaries lead to convolutional integrals and do not therefore destroy 
the O(n) memory demand of the method. The FE-CGFFT formulation relating to circular 
(and ogival) boundary enclosures is described in the University of Michigan report 
025921-1 1-T (see also [13]) and results based on its implementation are shown in Figures 
5 and 6. Fig. 5 shows bistatic scattering patterns for a coated circular cylinder with a 
conductor radius of 3A, 0.05A coating thickness and material properties Cj = 3-j5 and |i r = 
1.5 - j0.5. The agreement with the series solution result is excellent. In Fig. 6 a 
backscatter pattern is shown for a A/2 x 1A conducting ogive (see Fig. 3). In comparison 
with the moment method results, the agreement is again excellent. Additional results are 
given in Figure 7 for a missle-like shape scatterer. 


Pre- and Post-Processing A Jgorithms 

The availability of pre- and post-processing algorithms is crucial for the generation 
of the geometry and display of results in a graphical form. Generally, it is desirable that 
these tasks be done with a graphical user interface (GUI) and possibly in an X-window 
setting. Part of this year’s effort was therefore devoted to the development of such 
algorithms and/or interfaces for the more sophisticated commercial pre- and post- 
processing packages. 
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For the most part, there exist commercial geometry, mesh generation and post- 
processing packages which are highly interactive and graphical. Nevertheless, there is 
always a need for a suitable interface or data interpretor between the commercial packages 
and the solution algorithm described in the previous section. The specific package 
interfaced with the computational algorithm was SDRC IDEAS@ and the selection of this 
was based on its availability on the U-M Network, its versatility, graphical user interface, 
and capability to generate meshes for 2D and 3D structures. Furthermore, a new version of 
SDRC IDEAS@, to be released soon, will support X-windows. IDEAS was developed 
for mechanical design purposes, but its geometry and finite element mesh generation 
modules are particularly suited for our needs. The geometry is defined graphically using 
the area (for 2D) or solids (for 3D) modeling capability provided by the module Geomod. 
Alternatively, the user may choose to enter the geometry in terms of individual points, 
surve segments (for 2D) or surfaces (for 3D). Once the geometry is entered, mesh areas or 
regions are specified and either a mapped mesh or free mesh can be generated. Further, 
individual nodes and/or elements may be inputted as desired using the CREATE command. 
Once the mesh is generated, two files are created, one containing the nodes and their 
corresponding coordinates, and another specifying the nodes of each element. These files 
are then read by an interpretor which creates a new input file compatible with the format 
required by the computational modules. 

Examples of two-dimensional meshes generated with SDRC EDEAS@ were shown 
in Figures 1, 2b and 7, and these are in a form suitable for the FE-CGFFT analysis. Some 
three-dimensional meshes are also displayed in Figure 8 for an ogive and missile-like 
structures. As seen, the 3D meshes are terminated at a cylindrical surface, tightly enclosing 
the scatterer which is the intended enclosure for the proposed FE-CGFFT method. A brief 
manual for geometry and mesh generation using SDRC IDEAS@ is currently being 
prepared. 

SDRC IDEAS@ is a rather sophisticated package and its use is certainly preferable 
for 3D modeling and mesh generation. For 2D mesh generation, though, it is possible to 
construct a non- proprietory package without much effort, and which is also simpler 
without a serious sacrifice in versatility. Clearly, the primary reason for resorting to such 
an algorithm is to permit mesh generation at sites not having a license for SDRC IDEAS@. 
The specific geometry and mesh generation package developed for this purpose is based on 
the algorithm described in [14], Examples of free meshes generated by this package are 
displayed in Figures 2a and 3. The package is interactive/menu driven and can be readily 
used without much preparation. The mesh can be displayed in the Apollo screens or a 
postscript file may be generated for display on other workstations. At present, 
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visualization cannot be done in an X-window but this capability is planned for early next 
year. 

A variety of post-processing capabilities have also been exmployed for a graphical 
display of the output data. The output is either in the form of echowidth plots as a function 
of observation and/or incidence angle or in the form of gray level field maps. Color instead 
of gray level field maps can also be generated at those workstations which support this 
feature. Generally, all echowidth plots are generated and displayed using standard 
software, and each workstation provides its own selection. To generate and visualize the 
gray level and color field maps we employed the public domain package GRAFIC. An 
example of a gray level plot is shown in Figure 9. This is generated from a postscript file 
and can thus be displayed on other sponsor machines. 


3D Algorithm for Bodies of Revolution 

Before extending the presented 2D formulation to scattering by arbitrary 3D 
structures, it is instructive that we first consider its implementation for a restrictive class of 
3D bodies. In particular, during this year an algorithm was developed for scattering by 
inhomogeneous bodies of revolution. Because of the symmetry of this structure, it is only 
necessary to discretize it in a single plane slicing the structure as shown in Figure 10. A 
knowledge of the fields over this cross-section is then sufficient to generate the fields 
everywhere by employing a Fourier expansion in the azimuthal direction. Clearly, the 
discretization can be accomplished using a 2D mesh generation routine and this is the 
primary reasons for considering this class of structures. Also, the storage requirements are 
comparable to that of the 2D formulation although, as expected, the computational intensity 
is much greater. 

The mathematical details pertaining to the BOR formulation will be presented at 
sufficient detail in an upcoming technical report. Briefly, the method consists of the 
following steps 

1) A Fourier expansion is used to expand the fields in terms of those over a single 
cross-section of the BOR. 

2) The fields in the finite-element region are then formulated via the Coupled 
Azimuthal Potential (CAP) method as described in [15]. This results in a banded 
finite-element matrix in terms of the boundary fields. 
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3) The boundary fields are formulated via the usual Stratton-Chu equations which are 
then discretized via the boundary element method. As before, the boundary enclosure 
is chosen to yield convolutional integrals computed via the FFT. 

4) The finite-element and boundary-element systems are coupled via the boundary 
fields and solved via the CGFFT method maintaining an O(n) storage requirement, 
where n is the number of nodes over a single cross section of the BOR. 

Presently, a code has been written based on the proposed formulation and is in the 
final stages of the validation process. 


CONCLUSIONS 

o far, we have formulated and implemented the FE-CGFFT method for a variety 
of 2D structures and we are now in the process of completing its implementation for BOR 
structures. The method was proposed because of its versatility, accuracy and low memory 
demand in comparison with other methodologies, and all of these attributes have been 
demonstrated in the testing and validation process. It is therefore a promising method for 
general 3D implementations to be considered in the following year. 


TRANSITIONS 

The validation of the 3D BOR formulation is expected to be completed by early Fall 
1990. We will then begin the development and implementation of the formulation for 
arbitrary 3D structures. This implementation is expected to be much more involved than 
those considered earlier and the same if true for the geometry and mesh generation. It is 
therefore, likely that the proposed 3D implementation may not be completed by the end of 
the 3rd year. Also, because of the need to generate suitable pre-processing and post- 
processing algorithms additional man-hours are required during the third year of this effort. 
Most likely, a practical user-oriented validated and benchmarked code will not be available 
until the fourth year. As part of this effort it would also be desirable to design and develop 
a graphical user interface (GUI) compatible with the X-window platforms. The GUI is 
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particularly necessary for the 3D analysis package. Otherwise, the user would be faced 
with a long list of subprograms whose interfacing would likely be cumbersome. 

The proposed FE-CGFFT formulation employs an exact boundary condition at the 
termination of the me k h . This eliminates a need to extend the mesh far from the scatterer 
leading to a substantial savings in storage requirements. However, this storage reduction 
and solution accuracy is achieved at the expense of computational complexity and intensity. 
In many cases, though, where accuracy is not of primary concern, one could resort to the 
use of non-exact (i.e. absorbing boundary conditions), for terminating the mesh. This 
leads to completely sparse matricies which can be solved more efficiently using special 
purpose algorithms. In the future, it is therefore desirable to include this formulation as an 
option to the user. Also, a new class of boundary conditions are currently being 
investigated for terminating the mesh. 
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Figure 1 : E z and H z backscatter echowidth patterns for the illustrated 

perfectly conducting cylinder. 
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Figure 2(b): Example of a rectangular mesh enclosing a coated ogive. 


FILE: ogive Jeo_out 
CREATION DATE: 7 Aug 1990 
STRUCTURE: 0.5 x 1 X coated conductor 
ENCLOSURE: rectangle 

number of nodes: 1084 

number of elements: 1936 

nodes on pec boundary: 140 

nodes on obs boundary: 92 

nodes on unknowns (E-pol): 944 
nodes on unknowns (H-pol): 1084 
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Figure 2(c): Fl-polarization backscatter pattern for the coated ogive in 

Fig. 2(b). 
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Coated Conducting Cylinder R=3 X (H-pol) 
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Figure 5: E z and H z bistatic echowidth from a coated circular cylinder with 

a conductor radius of and coating thickness of 0.05X with 
material properites e r = 5 - j5, p r = 1 .5 - j0.5. 









FILE: missi!e_out 
CREATION DATE: 1 Aug 1990 
STRUCTURE: missile shape 
ENCLOSURE: rectangle 

number of nodes: 636 

number of elements: 908 

nodes on pec boundary: 172 

nodes on obs boundary: 1 92 

nodes on unknowns (E-pol): 464 
nodes on unknowns (H-pol): 636 
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Figure 7: Backscatter patterns for a missile-like perfectly conducting and 

coated cylinder. 


























Figure 10: Body of revolution 


surrounded with a rectangular mesh. 
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ABSTRACT 

The diffraction by a material discontinuity in a thick dielectric/ferrite layer 
is considered by modelling the layer as a distributed current sheet obey- 
ing generalized sheet transition conditions (GSTCs). The sheet currents 
are then formulated and solved via the standard dual integral equation ap- 
proach. This yields the diffracted field in terms of unknown constants which 
underscore the non- uniqueness of the GSTC current sheet representation. 
The constants are dependent on the geometry and properties of the discon- 
tinuity and are determined by enforcing field continuity across the material 
junction. This requires the field internal to the slab which are determined 
from the external ones via analytic continuity. Results are given which 
validate the solution and demonstrate the importance of the constants. 


OBJECTIVE 

This task involves the use of higher order boundary conditions to generate 
new solutions in diffraction theory. In particular, diffraction coefficients will 
be developed for dielectric/magnetic layers and metal-dielectric junctions 
which are often encountered on airborne vehicles as terminations of coatings 
and conformal antennas. Solutions for both polarizations will be developed 
for fairly thick junctions and versatile computer codes will be written and 
tested. Creeping wave diffraction coefficients will be also developed for 
multilayered coated cylinders. 

PROGRESS 
1 Introduction 

In scattering, layered materials are often modeled by equivalent sheets 
satisfying simple boundary/transition conditions. In particular, impene- 
trable layers are typically replaced by opaque sheets satisfying standard 
impedance boundary conditions (SIBCs) [1], whereas penetrable layers are 
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represented by transparent sheets obeying resistive or conductive type tran- 
s: ion conditions (STCs) [2], These simple boundary/transition conditions 
re ate the normal fields to their first normal derivatives through proportion- 

!f in the SIBC case and “resistivities/conductivities” 

in le STC case. With this modeling scheme, a discontinuity in layered ma- 
is represented by an equivalent sheet discontinuity, whose scattering; 
may be treated via function theoretic techniques such as the Wiener- Hopf 
method As is well known, however, these equivalent sheet representations 
are valid only for very thin or lossy layers and alternative simulations are 
therefore required to model discontinuities in low loss layers and/or layers 
of appreciable thickness. 

i . 0116 ,™ 1 ^ P r F w Ch iS t0 empl ° y S enera fized impedance boundary con- 
* 1 ° I ] 1S (^IBCs) [3] [4] or generalized sheet transition conditions (GSTCs) 

I5j [6J in place of the usual SIBCs and STCs. The GIBCs and GSTCs are 
respective generalizations of SIBCs and STCs and permit a more accurate 
representation of the fields at the surface of the coating or layer Unlike 
the SIBCs or STCs, GIBCs and GSTCs include second and possibly higher 
order derivatives of the field components on the equivalent sheet which are 
responsible for the higher accuracy of the conditions. The highest deriva- 
tive kept m the condition defines their order and generally the accuracy of 
the conditions is analogous to the order. As can be expected, thicker 'and 
multilayer coatings require higher order conditions for an accurate simu- 
lation and to date a plethora of GIBCs and GSTCs have been derived to 
model a variety of material coatings and layers [6] [7] [8] [9] [10]. 

GIBC/GSTC sheets are well suited for characterizing the diffraction 
by discontinuities in thick coatings or layers. In particular, they can be 
employed in conjunction with the Wiener-Hopf method or dual integral 
equation approach without much deviation from the procedure used in 
connection with the SIBC or STC conditions. However, the resulting so- 
lutions obtained in this manner are inherently non-unique [11] [12] This 
non-uniqueness cannot be removed with the usual application of the edge 
condition or the enforcement of reciprocity, which has been used in the past 
to generate a more physically appealing, if not a unique, solution 

L mqueness is an obvious requirement of the physical problem and unless 
resolved it would seriously undermine the usefulness of the conditions. In 
he case at hand, the non-uniqueness is manifested in the form of unknown 
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solution constants [12] and this simply points to the fact that additional 
conditions are required for their specification. 

In this paper we demonstrate that the GIBC/GSTC sheet characteriza- 
tion can yield a complete solution when supplemented with certain condi- 
tions at the sheet discontinuity which do not require apriori knowledge of 
the edge fields. As a vehicle in presenting this solution procedure we employ 
the dual integral equation method to consider the plane wave diffraction 
by a discontinuous distributed sheet (see Figure 1(b)). This very general 
model is capable of representing material half-planes, material junctions, 
and material discontinuities on grounded structures, such as those shown 
in Figure 2. In addition, a distributed sheet model typically renders the 
same degree of accuracy as the usual infinitely-thin sheet, but with a lower 
order condition. It is, therefore, of much practical interest. 

In the first part of the paper, the GSTC representation of the distributed 
sheet discontinuity is used to develop dual integral equations in terms of 
the unknown spectral functions proportional to the sheet currents. These 
equations are then solved in the standard manner to yield expressions for 
the spectral functions in terms of unknown com ants, and examples are 
presented where a proper choice for the constants demonstrates that they 
recover known solutions. This demonstrates the validity of the presented 
solution, but in general, the determination of the constants requires the 
enforcement of additional constraints demanding field continuity across the 
junction. The development of these conditions and their use in solving for 
the constants is also presented. 

2 Dual Integral Equation Formulation 

Consider a distributed sheet of thickness r illuminated by the plane wave 

p _ Jk(r cos* 0+! ,si„<M = f E * Polarization, 

inc \ Z 0 H Zttnc , H z polarization, 

as shown in Figure 1(a). The excitation (1) induces reflected and trans- 
mitted fields which are explicitly given by the properties of the distributed 
sheet. If this sheet models a symmetric slab, then an appropriate GSTC 
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representation is formally given by [10] 


(-£) i F+ - F -} + (~£) {* [^ + + *“]}- 0, 


-oo < X < oo 




k "22 


dx 2 


k 2 


){aj,[F + -F-]} = 0 , 

— oo < x < oo. 


( 2 ) 

in which F is the total field, F± = F(x,y = ± r /2), flrF* = J. F (x y 

= ±r/2), and dyF* = f y F(x, y) | v=±r/2 . Also, ^ (-M) are differential 
operators which operate on the field quantity in the curly brackets, and are 
mte polynomials in whose coefficients depend on the slab modeled 

by the distributed sheet. To maintain the generality of the solution, the U l 
operators are left in symbolic form and the reader is referred to [10] for their 
explicit representation in terms of the material constants and thickness of 
the ayers comprising the modeled slab. In general, the order of U\ x (i e 
the highest derivative present) is usually the same or one more than that 

° S Slmilarl >' the order is the same or one more than the order 

ot U 22 . Thus, we may define the orders of the GSTCs in (2) to be 


]\° dd _ maximum { order of U\ x (A 2 ) , 1 + order of U\ 2 (A 2 ) } 


(3) 


Nl ven - maximum { order of U\ x (A 2 ) , 1 + order of U\ 2 (A 2 ) } 

The reflected and transmitted fields may now be easily determined by 

TN I rvt / O \ 1 _ 1 J 


employing (2) to find 


F re f[ = R ie i k ( xcos, t>o-ysin<f> 0 ) 

F lran = T’jgjfctxcostfo+ysin^o) 


m which R x and T, are the reflection and transmission coefficients 
tively, and are given as 

gjfcr sin i j> Q 

R i = — [R\ wn + R\ dd J 

e jkr sin0 o 

T i = — g [R\ ven - R° dd \ . 


(4) 

( 5 ) 

respec- 


( 6 ) 

(7) 
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with 


sin (cos 2 <j>o) - U\x (cos 2 <i>o) 

sin (PM 22 ( cos2 4>o) + Wjx (cos 2 <j) 0 ) 

sin <f>JA\ 2 (cos 2 <f> 0 ) - U\y (cos 2 <t> 0 ) 
sin (fiMh (cos 2 4> 0 ) + Uh (cos 2 <f> 0 )' 

We remark that in (8) and (9), U\ } (cos 2 <f > Q ) now represent simple polynomial 
functions in cos 2 <f> 0 , since —dx 2 /k 2 = cos 2 <p a in view of the field expressions 
(4) and (5). 

Consider now the case where the right half of the distributed sheet 
in Figure 1(a) is replaced by another sheet of the same thickness, but of 
different properties, as illustrated in Figure 1(b). The GSTC representation 
of this modified sheet is 

( 10 ) 

for — 00 < x < 0 and 

( 11 ) 

w| '(-^){ Ft+F ~} + i w -(-^){^[ F+ - F "]} = 

for 0 < x < 00 , where the superscripts 1 and 2 distinguish the left- and 
right-hand sheets, respectively. Referring to our previous discussion, the 
orders of the right hand side GSTCs are given as 

N odd 

— max {order of U\ x {A 2 ) in A, 1 + order of U x2 {A 2 ) in A j 

( 12 ) 

AT 6 " — max {order of U 2l (A 2 ) in A, 1 + order of U 22 (A 2 ) in Aj . 


R? 


R°\ 


odd 
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The modified right hand side sheet induces a scattered field F s in the 
presence of the excitation (1), and the total field can be represented as 


F = 


Fine + Frefl + Fs V > T / 2 
Ftran + F» V < r /2 


(13) 


where F 3 is the unknown scattered field in the region \y\ > r/2 and can be 
expressed as [13] [14] 


F 3 (x,y) = J c 


F t 3c id ( COS OL ) + Pt even ( COS Of ) 

y 


e -jks'ma(\y\-T/2) e -jkxcosa^ a ^ 


(14) 


where C is a contour in the complex a plane, such that A = cos a runs from 
— oo to oo as shown in Figure 3. In this, the spectral functions P 0 dd (cos a) 
and Peven (cos a) are directly related to the Fourier transforms of the un- 
known equivalent currents 


Jodd = F+ - F~ (15) 

J even = F+ + F ~ , (16) 

via the relations 

/ OO fl\ 

P odd (\)e-> kxX -=== (17) 

-oo y 1 — \ z 

J even (x) = 2 /^Peve n (X)e-^ x ~L= . (IS) 

Substituting (1), (4), (5), (13) and (14) into the transition conditions (11) 
and ( 12), and introducing the transformation A = cos a (see Figure 3) yields 

Q° dd (A 2 ) P od d (A) e -i- " = 0 (19) 

/” Qt ven (A 2 ) P even (A) = 0, (20) 
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for x < 0 and 


/ 0 ° , . 

2 G° 2 dd A 2 )P odd (X) 

-oo x f 


jkx\ 


d A 


v / n r A2 

2 sin <j> 0 e ]kxX °e jkT ' 2 * m4> °Z odd ( A 2 ) 


<?f“(A 2 ) 


/~ 2C?r en (A 2 ) (A)e- 2fclA 7 =£ 


vt^a 2 


[2i: 



2 sin <t> 0 ei kxX °e jkT P* in4, °Z even (A 2 ) 

(22) 


Q e r n { a 2 ) 

where A 0 = cos d> 0 and 


er"M = tt'n (a 2 ) + (a 2 ) 

(23) 

Qr n (a 2 ) = u\ x (a 2 ) + vr=rvu' 22 (a 2 ) 

(24) 

Qf l (a 2 ) = U 2 XX (a 2 ) + yj\ - A 2 Z4 (a 2 ) 

(25) 

£f' en (a 2 ) = U\ x (a 2 ) + JTxy>ui 2 (a 2 ) 

(26) 

= 1 


(27) 

7 f\ 2N ) - 1 

even l / ' 0 J — 

K (A|) (A2) - Wi (A J)w|, (A|)] 

(28) 


Equations (19) with (21) and (20) with (22) form two uncoupled sets 
of integral equations, sufficient to yield a solution for the unknown spectra 
Podd{\) and P e t , en (A). Clearly, because of the similarity between the two 
sets of equations, once a solution for P odd ( A) is found, the corresponding 
one for P even (A) follows by inspection. 
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3 Solution of the Diffracted Field 


Upon a solution of the dual integral equations (19) and (21) we obtain 

p jkr / 2 sin0 o 


7 S 1 I 10 

P -' A > - h—x 


- A 2 


A 0 GfJ (A) gfj (A 0 ) G° 2 f (A) gf + d (A 0 ) 


J odd 




Eodd (A) 
E odd ( A 0 ) 


N, odd 1 N 0 dd 1 

+ °mn 


(A + A 0 ) m (AA o ) n 


m=l 


n=0 


( 29 ) 


where we have assumed that J 0 dd{x) ~ x Sodd as x — ► 0 with 0 < S 0 dd < 1- In 
this, N od d = int {l/'2(N^ dd + N^ dd +1)}, and a mn are arbitrary constants 
as yet undetermined, and correspond to the coefficients of the polynomial 
resulting from the application of Liouville’s theorem. The chosen symmetric 
form of this polynomial is not unique but will be found most useful later in 
constructing a reciprocal form for P 0 dd{ A). Also, Gt+(A) are Wiener-Hopf 
split function regular in the upper (+) or lower (-) half of the A-plane and 
satisfy the relation 


Gi(A 2 ) = G 1+ (A)Gi_(A) 


(30) 


(see Appendix). Similarly, <7 2 ±(A) are the corresponding split functions 
associated with G 2 (A 2 ). Finally, E 0 dd( A) is some entire function behaving 
than \\\{ N l<ii +N ldd)l 2 ~ s o<id anc j can take an y 0 f the forms 

E 0 dd ( AA 0 ) or 


no worse 


E 0 dd( X ) — { 


Z od d( X ) or 

ZtddW 


(31) 


where Z ± ( A) are again upper and lower functions satisfying the relation 

Z( A 2 ) - Z + (A)Z“(A) (32) 

Following a similar procedure we obtain P even (A) as 

gjkr /2 sin 4> 0 


j sin <j> 0 yA - A 2 
(A) " 2?r 


A + A 0 Q\l en (A) Qf™ ( A 0 ) G e 2 v + en ( A ) Ql v + en ( A 0 ) 




Z 7 1 ( \ \ A Nevcn — 1 — m 

r , + E E <>„„(A + A X(Uo) n 

-Jeven \ 'V) m=l 


n= 0 


(33) 
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with E even (A), N even and b mn being the counterparts of E odd ( A), N od(i and 
a mn , respectively. Taking into account the choices (31), we may substitute 
(29) and (33) into (14) and subsequently perform a steepest descent path 
evaluation to obtain for p — ► oo (all surface wave contributions are neglected 
in this evaluation) 


e~ jkp 

F (p,<f>) ~ [D od d (<?b 4>o) + D even ( <f> , <j> 0 )\ - -y ■ - f 34) 

yjkp/2n 

where (p, <£) are the usual cylindrical coordinates and D odd (</>, 0 o )+Asven (0, <^ 0 ) 
is the far zone diffraction coefficient symmetric with respect to <b and <p Q . 
We have 


D odd ( 0 , ) 


e J7r ^ 4 sin <f> Q sin <f> 

2n cos <£> + cos 

e jArr/2(sin<£ o +| s in0|) 


(cos 0) C/f" (cos (f> 0 ) (cos 0) (cos <j> 0 ) 

' ~ Nodd - 1 N odd -l -m 

Z odd (COS <£>, COS <j4 0 ) -f" E E 

m = l n=0 


D even ( (f), <f> 0 ) — 


• (cos <p -f COS <f>o) m (cos <j> cos <f>o) n 

e -j>/4 s i n <^ o | s i n <*>| 


(35) 


27T COS (f> + COS 4> 0 

e j/cr/2(sin0 o +| sin0|) 


£ 1 - n (cos <f>) g\ v * n (cos <p 0 ) g%l en (cos <£) 0f£ en (cos 4>o) 

^ — 1 N cv«n 1 — 771 

% even (cos <j), cos <f > 0 ) T E E ^771 Tl 

m = l n=0 


(cos +COS ^o)" 1 (cos cos <f> 0 ) n 


(36) 
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( 37 ) 


in which the functions Z 0 dd,even (cos cos <f> 0 ) are given by (see (31)) 

i Zodd ( - cos </> cos <j>o ) or 

Z~dd ( cos <£) Z odd ( cos ^o) or 

Ztdd ( cos 4>) Z odd ( cos 4>o) 

{ Z ev en ( - COS <j> COS <f> 0 ) Or 

Z 7ven ( COS <j ) ) Z ~ en (cOS <f> 0 ) OT ( 38 ) 

Z Ln ( COS <t>) Z Ln ( COS <f>o) 

Because the above three choices for Z 0 dd and Z even differ only by terms of 
the form (cos (j> + cos <f> 0 ) m (cos <j > cos <j> 0 ) n , it is immaterial wdiich of them we 
choose, although one of the choices may likely lead to a more compact rep- 
resentation. Nevertheless, regardless of the choice of Z Q dd and Z even , one is 
still faced with the determination of the unknown constants a mn and b mn in 
(35) and (36), repectively. These are a manifestation of the non-uniqueness 
of the finite-order GSTC sheet model employed herein, and their explicit 
determination requires the introduction of additional constraints pertain- 
ing to the physics of the problem. Before we consider their determination 
for the general case, we first look at a specific example, that of diffraction 
by a thin single layer junction. 


4 Diffraction by Thin Single Layer Discon- 
tinuous Slabs 

The diffraction coefficient given by (35) and (36) is very general and can 
model a wide variety of geometries. To check its validity, display its ver- 
satility, and assess the relative importance of the unknown constants, we 
consider the thin material-to-material junction of thickness 2 w as shown 
in Figure 4. The slab will be modelled by a sheet of thickness 2 (w — w s ) 
and with a proper choice of the material parameters this geometry can re- 
duce to junctions whose diffracted field is available, thus, permitting some 
validation of our solution. 

If the left hand side of the slab, in addition to being thin, is also asso- 
ciated with low index of refraction, it may be modeled by a low contrast 
GSTC sheet. Thus, an O(u) 1 , w]) approximation with terms of O(uqw) ne- 
glected is sufficient for the representation of the operators or polynomials 
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1 A \ . In particular, we have 


U\ x (- dx 2 /k 2 ) = 1 

U{ 2 (-9.r 2 /^' 2 ) = jk(uiw-w a ) 

U l n (-dx 2 /k 2 ) = jk - , 

(~dx 2 /k 2 ) = 1 


- w *) + i U ” Ws ) 9x2 


where e x and /./ i are the relative permittivity and permeability of the left 
hand slab, respectively, and 

( fi i, polarization 
1 \ e x H z polarization 

Also, when w s = ic, these are simply the transition conditions derived first 
by Weinstein [5] and later by Senior and Volakis [6]. The corresponding 
polynomials to be employed in (23) - (28) are given by 

U{\ ( — cos (f) cos <f > Q ) = 1 

U\ 2 ( — cos <$> cos (j> 0 ) = jk ( u i w — w s ) 

UL ( — cos <p cos (f > 0 ) = j k ( — — w s ) + j k ( w s ) cos <f> cos d> 0 

\ Ul J \U\ ) 

U 22 ( — cos <f> cos <f> 0 ) = 1 (41 ) 

Incorporating these into (35) and (36) and setting 

z odd (cos 0, cos <j> 0 ) = Z 0 rf t /(-cos^cos^ 0 ) (42) 

Zeven (COS (j)^ COS (f ) 0 ) — Z GV€n ( COS (f) COS ^> 0 ) (43 ) 


Z^odd ^o) — 


J7r / 4 sin0 o sin^ 

< 

2tt cos^ + cos < f > 0 


p jA:r/2(sin <£ 0 + | sin<£|) 


U \ 2 ( — cos (j) cos <j> Q ) — j k (u\w — w s ) U xl (— cos <f> cos <f> 0 ) 
M- (cos <p\ 7 odd ' x ) M- (cos <j ) 0 ; 7 odd fl ) £2+* ( cos <£) ( cos 0 o) 
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& even ^o) — 


sin 4>p\ sin 0} c ;fer/2(sin 0 O + 1 sin0|) 
2;r cos0+cos0 o 


[a 1 ~r U2 ms <Po J t - 1 22 V y 

a3 [rim=l ( cos ^ i7m * ” ’ 1 )M-(c08(^o;7m' 


Vj y 2 + 


+ 


6 1 o (cos 0+cos 0 O ) 




l\»r { even,l\ 

y\f- {COS VoWm 


)^ en {cos0)^|; en (cos0o) 


where the split fui^ f ‘on Af_ (cos <£ 57 ) is given in the Appendix, 


and 


OC\ 

= jk 

(we i/zi 

1 W i 

\ U\ 

Ot 2 

= jk 

(“»*) 

<*3 

jk 

-(ei/ii-1) 


u 1 


7 


o<id,l 


-3 

k (uiw — w s ) 


( 45 ) 


(46) 


even. 1 

7l,2 


til ± \Ju\ + ik 2 w (ei^i — 1 ) (u> — U\W S ) 
2 jk(w - w s ui ) 


(47) 


with 7 0<w OT ei ’ en are associated with possible surface wave poles. To complete 
the definition of ( 44 ) and ( 45 ), the functions associated with the right hand 
side properties of the slab (i.e. those functions with the superscript 2 ) must 
be specified and Tables 1 and 2 provide explicit expressions for the functions 
Uf- ( — cos <f> cos <j> 0 ), g& d (cos4>)Q° 2 d + d (cos £,) and Q™ n {cos<f>) Q% en (cos <j> 0 ) 
terms. By edge condition considerations, all of the constants a mn and b mn 
have been set to zero except bio appearing in the definition of D even , which 
is non-zero unless the right hand side slab is a PEC/PMC under an E z / H z 
excitation (see Table 2). 

By invoking image theory the diffraction coefficient for the 
grounded metal- dielectric join, shown in Figure 5 is given by 

D rs { cos <j), cos <t>o) = 2 D euen (cos </>, cos 4>o ) (48) 
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The GSTC or GIBC model for this structure cannot dsicriminate whether 
the stub at the junction is a perfect electric conductor (PEC) or perfect 
magnetic conductor (PMC). This information can only be carried by the 
constant i 10 and its determination must somehow involve the properties 
of the junction across its thickness as discussed in a subsequent section. 
However, since the diffraction coefficient for the junction in Figure 5 is 
already available [15], 6i 0 can be identified. Upon setting w s = 0, we find 



This comparison clearly demonstrates the importance of the constant 6 10 
and by referring to Figure 6 we observe that it plays a major role in the 
computation of the diffracted field. 

5 Modal Decomposition of the Symmetric 
Slab Fields 

A general approach for determining the solution constants is to enforce 
tangential field continuity accross the junction. This, of course, demands 
a knowledge of the fields internal to the discontinuous slab, which are not 
readily available when a GSTC simulation is employed. The Weiner- Hopf 
(or dual integral equations) solution in conjunction with the GSTC pro- 
vides only the fields external to the slab, and this section deals with the 
determination of the internal fields from the external ones. 

A modal representation of the internal field is first proposed comprised 
of discrete and continuous spectral components. This representation is 
compatible with that given by Shevchenko [16] whose eigenfunctions are 
chosen to satisfy field continuity across all layer interfaces including the 
air-dielectric interface. Consequently, the representation is valid inside and 
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outside the dielectric once the coefficients of the modal representation are 
determined. This is accomplished by recasting the Weiner-Hopf or dual 
integral equation solution given earlier (see (13), (14), (29) and (33)) in a 
form compatible with the proposed modal representation, thus permitting 
the identification of the modal or eigenfunction coefficients. These will, of 
course, be in terms of the unknown constants appearing in the Weiner-Hopf 
solution and the enforcement of field continuity acoross the junction leads 
to a linear system of equations to be solved for the constants as described 
in the next section. 

For the symmetric slab in Figure 2, the total field may be decomposed 
into its odd and even components. Specifically we write 

f F l ' odd (x, y) + F 1,even (x, y) x < 0 
F(x,y) = { (52) 

{ F 2 ’ odd (x, y) + F 2 ' even (x, y) x > 0 

where F odd {x,y) = —F odd (x,—y) and F even (x,y) = F even (x, —y). Follow- 
ing [16], the odd and even fields interior and exterior to the slab may be 
expanded into discrete and continuous eigenmodes as 

N go s 9 \ 

41, odd ^1 , odd / ^1 ,goy \ e -jkxy„\ 9 ° 

m- 1 ' ‘ 

AT 1, Odd 

1 3W 1 

+ £ B l m odd ^ m odd (y)e-^ x - 

m= 1 

+ r C hodd (/?) V hodd (a 2 , y) e~ 3kxX d(3 (53) 


~il .odd 


(avi/) = 


F 


1 ,even 


(*,y) 


-jkx\ 


m= 1 V 7 

Krl,even 
1 v aw 

+ E B^ even ^ ven (y)e~ 3kxX 

m = 1 

+ C 1 ' even (/?) *F 1,even (A 2 , y) e~ 3kxX d/3 


\ 1 ,cven 


( 54 ) 


.odd 


(*>y) 


Ngo j \ 

= $;/!«'♦« ()?•) ,y) 

m—1 


-jkxX 


2, go 
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*r2,odd 
1 y SW 


,2, odd 


+ E B 2 m 0dd e m » dd {y)e-> kxX " 

m= 1 

roo . 

+ J C 2 ’ odd (/?) <H 2 ' odd (a 2 , y) e~ jkxX d/3 


( 55 ) 


E ven ^ 2 ’ euen ((A^r) 2 , y ) 

m = l ' ' 

»r2 ,«ven 

™4tl/ 

+ E rT ' (y) e~i kxX ™ V * n 

m — 1 

zoo 

+ / C 2 ' even (j3) * 2 ' even (a 2 , y) e~* kxX d8 (56) 

where Im{\° dd,even } < 0 and A = \/l — d 2 , with the branch of the square 
root chosen so that Im {^/ 1 - /? 2 } < 0. In (53)- (56), 'P 0(RcveTl are referred 
to as the cross section functions corresponding to the continuous modal 
fields whereas $^’ et;en are the corresponding cross section functions for the 
discrete modal fields associated with the surface waves. The cross section 
function associated with the geometrical optics fields is also \fy odd ’ even eval- 
uated at A = A^\ where A^is a parameter to be determined later. As can 
be observed from (53) - (56), the cross section functions specify the field 
behavior in the plane normal to the slab, and hence all information per- 
taining to the fields interior to the slab are embedded into these functions. 
They will be chosen to satisfy the orthogonality relations (where u.(y) is 
/i(j/) or e(y) for E z or H s polarization, respectively) 


r<2,et ;en 


(*,y) = 



tt(A 2 ,y)fr(A 2 ,y) 

u(y) 


dy 



*m(y)*(X Vy) 

u(y) 


dy 


0 for A ^ A 

0 


(57) 

(58) 


and thus each discrete eigenmode (y) e~ jkxX and each continuous eigen- 
mode ^ (A 2 , y) e~j klX must satisfy the wave equation. Additional details 
pertaining to the cross section functions are given in [16]. 
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Exterior Cross Section Functions 


To compute the cross section functions in the exterior slab region |y| > 
r/ 2 , we recall that in accordance with the slab simulation based on the 
generalized sheet transition conditions (GSTCs), the external fields satisfy 
the conditions ( 10 ) and (11). Because of the orthogonality relations (57) 
and (58), each of the cross section functions (A 2 , y) and (y) must then 
satisfy their respective odd or even GSTC. In view of this we set 



+Wf 2 (A 2 )cos k(\y\ - r/ 2 ) Cl - A 2 ]} 


= y {Uu (A 2 ) ?i 2 (l, 1, |y| - r/2, A 2 ) 

+U P n (A 2 ) 922(1, 1 , |y| - r/2, A 2 )} (59) 



+U% 2 (a 2 ) cos [A: (|j/| - r/2) - A 2 ] } 


= {^21 (A 2 ) 912 ( 1 , 1 , |y| — r/2, A 2 ) 

+W 2 P 2 (A 2 ) 922(1, 1, |y| - r/2, A 2 )} (60) 

where 9 ^ represents the infinite order form of the 9 ,j layer operators given 
in [10]. Once each of the modes comprising (53) through (56) is substituted 
into ( 10 ) or ( 11 ), the differentiation implied by —dx 2 /k 2 reduces to a mul- 
tiplication by A 2 and the above T p ’ 0<w and 'I ,p,et ' en are then readily shown 
to satisfy the associated GSTC. It can also be shown that these satisfy the 
orthogonality conditions (57) and (58). 

A customary representation for the surface wave cross section functions 


is 


K odd (y ) = 

*r en (y) = 


\y\ e -im-r/2)^-(x^ d ? . | y) > r/2 
y 

e -MM-T/7)yJi-{Xr~) 2 . |y| > r/2 


(61) 

(62) 
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where \p^ even ' Q(i(i must now be chosen so that they satisfy their associated 
GSTC. By substituting (61) and (62) into (10) and (11), we find that 
^odd,eve n mus t satisfy the polynomial equations 




Kn 


odd 


+ Uf 


H \ rn 




,odd 


') - 


(63) 


^1 - (A ([A^ even ] 2 ) +U p n ([A^ l ' en ] 2 ) = 0 (64) 

and can be also identified as the poles of the slab plane wa v e reflection 
coefficient. We further note that 

l2 


fyPiOdd l hp 


$ 


,odd 


p.even 


Kodd 


,y) 


(y) = 


|y| > t/2 

(y) = 

^p,euen ^p,et-enj2 ^ ^ 

\y\ > T /' 2 

Ul, ([A5i*"f ) 


(65) 


( 66 ) 


implying that for a multilayer slab the cross section functions associated 
with the discrete and continuous eigenmodes are of the same generic form 
given by (59) and (60). 


Interior Cross Section Functions 


We consider now the determination of the cross section functions for the 
region interior to the slab (i.e. in the region |y| < r/ 2). For simplicity let 
us first assume a single layer slab of thickness t — T\, whose upper face is 
located at y = — Tj/ 2. In accordance with the preceeding, the cross section 
functions associated with the external fields are given by 

^■ M (k 2 ,y) = ld{?„(.<;.Kf,rf,A 2 ),- 12 (l,l,|y|-r 1 /2.A 2 ) 

+ yi 2 A 2 ) y 22 (l, 1, |y| - Ti/ 2,A 2 )} ; |y| > Tj/2 (67) 

* p ’ even (A 2 ,y) = {^(^^r^A^^iaJyl-n^A 2 ) 

+<? 22 («!,«!, rf, A 2 ) 922 ( 1 , 1 , |y| - ti/ 2,A 2 )}; |y| > n/2 (68) 
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obtained by setting U^-{\ 2 ) = («i, rf, A 2 ) in (59) and (60). These 

are orthogonal functions and each must, therefore, satisfy the continuity 
conditions 


yp,°dd ^2 r - 

) = <H p ' odd (\ V+) 

(69) 

— dy<a pMd (a 2 , t," 
u 1 v 

) = dyW' odd (A 2 ,r+) 

(70) 

<S> p ' even (A 2 , t~ 

) = <f' p ’ even (A 2 ,t+) 

(71) 

—dy^ p ' even (A a ,rf 

) = dy<H p ’ even 

) (72) 

with similar conditions on $£f i ' even (y). It is now straightforward to deduce 
that possible cross section functions satisfying (69) - (72) are of the form 

W odd (\ 2 ,y) = 

-- — yi2(«?,«?,|y|,A 2 ) 

y 

(73) 

V p ' even (A 2 , y) = 

-- q 22 (u p 1 ,K P 1 ,\y\,\' 2 ) 

(74) 


for \y\ < r/2. Also, in view of (69) - (72), the cross section functions for 
the surface wave modes remain as given in (65) - (66), provided (73) and 
(74) are used in place of \fj° dd ’ €ven m 

For the general case of a multilayer slab, it is necessary that each of 
the internal cross sections functions satisfy the continuity conditions at all 
layer interfaces comprising the slab. In this case we find that 


qjP>odd 



M 


Wn(A 3 )$ ia (l,l,M -r/ 2 , A 2 ) 

+U\2 (A 2 ) <722(1, 1 , \y\ ~ r/ 2 , A 2 ); 

|y| > r /2 

'Pu (A 2 )?i2(uf,K;f, \y\ - y/-i, A 2 ) 
+V['% (\ 2 )q 22 (tf,rf, |y| - yi-i 1 A 2 ); 

yi > \y\ > yi-i 

9i2(u?,K?, |y|, A 2 ); \y\ < y x 


(75) 
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^ 2 i(A 2 )?i 2 (l,l,|y|-r/ 2 ,A 2 ) 

+U% 2 (A 2 ) <122(1, 1 , |y| — t j 2 , A 2 ); 

|y| > r /2 

* p - even (A 2 ,y) = V21 (A 2 )yi2(«f,Kf, |y| - y/-i, A 2 ) ( 76 ) 

+^ > 22 P (A 2 ) y22( u f, «f, jy | — y/-i, A 2 ); 

yi > \y\ > yi-i 

. Q22(u p u kI, |y|,A 2 ); |y| < yi 



When these are used in (53) - (56) in conjunction with (65) and (66) we 
have a complete field representation for all x. 


6 Recasting of the Dual Integral Equation 
Solution for a Material Junction 

The expressions (53) - (56) can be used to represent the fields interior and 
exterior to the slab. It remains to find the coefficients of these expansions 
and to do this we must first rewrite F(x,y) in a form compatible with (53) 
- (56). That is, we need to identify from (13) and (14) the discrete and 
continuous spectral components. The discrete portion of the spectrum is, 
of course, comprised of the geometrical optics and the surface wave fields. 
These can be identified by detouring the integration path in (14) as shown 
in Figure 7. In particular, for x < 0 the integration path may be deformed 
to one over the branch cut in the upper half of the A plane, capturing 
any surface wave poles attributed to the zeros of (A) and {?[l en (A). 
Similarly, for x > 0, the integration path may be deformed to one over 
the branch cut in the lower half of the A plane causing the capture of the 
geometrical optics pole at A = — A 0 in addition to any surface wave poles 
attributed to the zeros of (A) and Q |^ cn (A). 
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Through the above deformation of the integration paths in (14) we 
obtain 


F° d \x,y) 


F even (x,y) 


F l g f dd (x, y) + F^ odd (x, y) + F d \f/ (x, y) x < 0 
Fgo° dd (*,») + F]f d (x, y) + Flff (x, y) x > 0 

F i g r n (*, y) + F s \-r n (x, y) + F^ v ; n (x, y ) x < o 

F g f ven (x, y) + F 2 ^ en (x, y) + Ftf™ (x, y) x > 0 


where the components F go , F sw , Fdiff denote the geometrical optics, surface 
wave, and branch cut (or diffraction) contributions to the total fields. 
After some manipulation we find 


fU ,odd 
30 


(x,y) 


A\' odd (A 0 ) V l ' odd (Xl y) e ]kxco * 4> ' 


(SO) 


,odd 

90 


x,y) 


Al' odd (A 0 ) ^ 2 ’ 0dd 



gjkx cos <£ 0 


(SI) 


F g f ven (x, y) == .4i' e,Jen (A c ) (A*, y) e jfc * cos *° 

(82) 

Fgo Ven (x,y) = Al' even (\ 0 )V 2 ’ even (\ 2 0 ,y)F k * COS(t> ° 

(S3) 

where the *4 expansion coefficients are identified as 


ii 

o 

-e 

~xi 

O 

sin (j) 0 eA T / 2s ™*° 

(84) 

w dd m 

A\' odd { A 0 ) = 

sin <j )o e jkT Fs“'*° 

(85) 

Gi dd ( AS) 

A}’ even (A 0 ) = 

sin ^ oe iW2 

(86) 

er n ( A o) 

j4i’ et,en ( A 0 ) = 

sin <t> 0 e jkT Fs™*o 

(87) 

Qeven (A 2) 

For the B coefficients, we have 



Ap-odd 

r 


Flf\x,y) = £ B,‘^(A„) 

/— 1 

^3 

"0 

O 

o 

~o 

1 
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Nodd-l N odd -\-m 


+ E E (A, + A„) m (a;-"a „)" 


m = l n=0 


,T # l,odei f \ xl.orfrfl 2 \ 

• * A / ■ i v e 


—jkx Aj ’ 


F'S' n (z,y) = E B l (A.) z„«(- A„A, ,,m ") 
/=i L 

Neven 1 ^cven — 1 — HI 


+ E E i m »(A; < “ n + A„)"(A; t ”'"A 0 )' 


m=l n=0 

,l,e„e„ /^f A l^enJ 2 \ c _jfc, A, 1 '*'" 


5/' 0(W (A c 


— sin <£ 0 

\ l r odd + a c 


ac^f.v) 

ax 


^jkr /2 sin 0 O 


, (Ao)$ff (xj- odd ) GS+(X 0 ) 


— sin <£ 0 

A l^et'en . \ 

/ i 


c jkr /2 sin 


v«n/ \ il 7 i \ A. ^ 7 

y l- v ; f?euen / \ \ /^eren f \ l,et/en\ ^ e z;en / \ \ 

aA J A = A 1.«*« ■ l - ^ °^ 2 + ^ 


The expressions for F?£ dd (x,y) and T; iu e, ' en (.r, y) parallel those in (86) and 

(87). 

To obtain the C coefficient we express F^fff as 

pi .odd, x y f°° j ~N l - A o 


r,i,odrf/ \ y f°° j -pyjL-A; 

d,ff X ' V |l/| /-oo 2?r _ £2 - /?2 + A 0 ) 

e jfcr/2( v /TTf+/j) e _jfc|j,| / J e _jfc Xv /i_ / 32 

£,°l d (v / w 1 + «Jj) 


7 (Ao) 
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Nodd-l N odd — 1 — 


✓ \ odd- 1 ; * od d ~ * “ 

■ - p) + E £ 

V ' m=l n=0 

•a mn (a o + \fl^Y (Xoy/l-^Y 


dp 


(92) 


where the branch of the square root is chosen so that L >/l — /i 2 ) > 0 and 
<5 is a vanishingly small positive number. By splitting the integral into its 
positive and negative portions, and employing some identities (92) becomes 

F\t,f(x,v) = 1-/J 2 ,A„) 

. $1 .Odd ^ ^ ^ e -,kX^P d/3 ( 93 ) 

with the expansion coefficient C 1 ' odd ({3 ) given by 


q\ , odd, fix _ l_ ^ 

(/ ' j -,yr^(yr^ + A 0 ) 



e }kTt2yJl- \%g 

odd | 

1+ 1 


1 

(~>odd | 
^2+ 1 

(v/l - P‘) S&MSitH K) {[M, 1 , (1 - /? 2 )] 

2 -P\u{ 2 (i-/? 2 )] 2 } 


Similar expressions for F%ff(x, y ), F] i ’jy n (x, y ) and Fj^ en can be obtained 
in a parallel manner leading to the identification of the remaining C coef- 
ficient. 


7 Determination of the Constants 

To determine the constants a mn and 6 mn , we may now enforce the tangential 
field continuity conditions 

F (x = 0 - ,y) = F (x = 0 + ,y) ; |y| < r/2 (95) 

^-arF (x, y) r=0 - = ;^y&zF (*> y).=o + ; lyl < T / 2 (96) 
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with 


“1,2 (y) = 


Mi, 2 (y) E z - pol 
^l, 2 (y) ^,-pol 


(97) 


and the subscripts 1 and 2 denoting quantities attributed to the left and 
right side of the slab. Substituting (79) - (79) into (95) and (96), we obtain 



77»2 ,odd 
go 


0-,«) + Fir" (* = 0-, y) + Fffi (x = 0-, y) 

(x = 0 + , y) + F%« (x = 0 + ,y) + Fl;ff (x = 0 + , y) 


(98) 


-i-rdx (*,») + Fit" (*, y) + Fj.vT (x, y)] x=0 . 


«2 (y ) 


<9x 


Fgo° dd i x ' y) + Fir“(*,y) + Flff (x,y)j x=o4 


(99) 


1 ,euen 
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T 72 .even 


0-,y) + Fir* (* = 0-,y) + Fi;7" (x = 0-,y) 

(x = 0 + , y) + (x = 0 + , y) + F%f{x = 0 + ,y) 


( 100 ) 


^9x [Fi— (X, y) + Fir™ (x. y) + eb7 ” (x,y)] i=0 _ 

= [F^“" (X, y) + Fir”*" (X, y) + F^7" (x, y)] i=0+ (101 ) 


to be solved for all a mn and b mn . In particular, for an odd GSTC of 0(N° dd ) 
to the left and of 0(N% dd ) to the right of the discontinuity, the number of 
a mn to be determined is equal to 


N„ 



N odd + N odd is even 
Nld d + Nldd is odd 


(102) 


To determine all a constants, (99) and/or (100) must then be enforced or 
sampled at a minimum of N a points accross |y| < r/2 and 0 < <f> 0 < tt. 
Similarly for an even GSTC of 0(N* ven ) to the left and of 0(N% ven ) to the 
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right of the discontinuity, 


N b = 


/V ( 1\J 1 

1Y even l^ v even 

2 

cvcn~^^cv<sn){^even'^^even~^ S } > iyl » 7^2 ’ 

^ a *” 7 ^euen ' ** even 1 


is even 


( ) “1 . 
8 


N l + N 2 

A piien. 1 ^ (°7n 


(103) 


is 


odd 


and thus, the b con «t ants can be determined by enforcing (101) and/or 
(101) at a minim, a of N b points. 

Substituting for the fields in (99) and (100) as given in the previous 
section, we obtain the equations 


V^(K,V) 

N a 

= ^2 a v Z°p d (m (p) , n (p) , A 0 , y) 

P= 1 

(104) 

V£¥( K,,») 

iV 

= J2 a pZdi d F (m(p) ,n(p) ,\ a ,y) 

p=l 

(105) 

where a p = a m ( p )„( p ) with 




(n + m — 1) (m + n) 

P = g +m 


(106) 


m{p) 


n{p) 


P 


-| /nt { ^±ybl } 


• Int 


= Int 


\/l + 8(p — 1) + 1 

2 

v /l + 8(p-l) + l 


m(p) 


(107) 

(108) 


which are in accordance with the ordering of the a mn constants as the order 
of the GSTC is increased (see Figure 8). The functions V F dd , V£ d p, Z F dd and 
Zg dd F are readily determined from the previous analysis and are not quoted 
here. 

Equations similar to (104) - (105) can be obtained for the b constants 
in a parallel manner. 


36 



8 Validation of the Solution 


The validity of the derived angular spectra and diffraction coefficients was 
already performed to a limited degree in Section 4 of the paper. What 
remains, therefore, is a validation of the procedure for computing the con- 
stants a mn and b mn which amounts to solving a small matrix. The valida- 
tion was done by comparison with processed data from a numerical model 
which consisted of a finite length slab having the prescribed discontinuity 
or junction at its center. First, the transient response of this finite slab 
was generated from bandlimited frequency domain data. The contribution 
from the material junction was then obtained by time gating the transient 
response. Numerically derived data from this procedure were found in good 
agreement with the presented analytical solution. An example is shown in 
Figure 9 corresponding to a thick (0.2 freespace wavelengths) material half 
plane. The numerical and analytical data are clearly in good agreement, 
and it is again demonstrated that the constants play a major role in the 
solution. 


9 Other Applications of the GIBC/GSTC 

It was shown above that the GIBC can effectively model thick planar layers 
of material. However, corresponding GIBC can also be derived for curved 
coated surfaces (see Figure 10), and the improved accuracy of these is 
particularly evident when surface wave effects are dominant. For surfaces 
having relatively large radii of curvature these can be easily derived from 
those of the planar surface with x and z replaced by the local tangential 
variables and y by the normal one. With a second order GIBC derived in 
this manner, the Mie series and GTD solutions have been found [20], [21] 
and compared with the exact modal series solution for a coated cylinder. 
As illustrated in Figure 11, the field given by the Mie series based on the 
GIBC is in excellent agreement with the exact result even at points close to 
the surface in the shadow region where a finite order boundary condition is 
inadequate. In contrast, data based on the standard impedance boundary 
condition (SIBC) are substantially inaccurate. 

Higher order boundary conditions have advantages in numerical treat- 
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ments as well. When used to simulate a coating, a GIBC eliminates the 
need to sample inside the dielectric, and this is important when storage is 
limited. In addition, it may be possible to use a GIBC to transfer a bound- 
ary condition to a plane, thereby producing a boundary integral equation of 
convolution type. In conjunction with an FFT, the equation can be solved 
iteratively to reduce the storage requirement to 0(n) where n is the num- 
ber of unknowns. As an example, for the three dimensional problem of a 
cavity in a coated ground plane, a GIBC provides a simple modal as well as 
a reduction in memory. If the coat^g is lossy or tapered in thickness, the 
non-uniqueness due to the terminations is avoided, and the same is true for 
cavities whose depth tapers to zero. Nevertheless, caution must be exer- 
cised when solving the integral equation numerically. The GIBC results in 
higher order derivatives applied to the Green’s function, and even if some 
can be transferred to the current, the increased singularity of the Green’s 
function makes discretization more difficult. In spite of this, integral equa- 
tion methods using GIBCs of up to the third oder have been successfully 
implemented [22]. 
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Appendix: Multiplicative Split Functions 

In this appendix we consider the splitting of 

Q (a 2 ) = U A (a 2 ) + vd - \ 2 U b (a 2 ) (109) 

as a product of two functions, one of which is free of poles, zeros and 
branch cuts in the upper half of the A plane and the other having the same 
properties in the lower half of the A plane. That is, we seek to write Q ( A 2 ) 
in the form 

Q (A 2 ) = a+(A)0_(A) (110) 

where the superscript + and - indicate an upper or lower function, respec- 
tively. Noting that 


U A (A 2 ) 

4 

= 

n — Q 

(111 

Ub (A 2 ) 

N B 

= EM 1 -* 2 ]" 

n— 0 

(112 


with N,\ — Nb or N A — Nb + 1, we may rewrite Q (A 2 ) as 

N s 

5 (V) = £ S n [%/i^v]" (113) 

n— 0 

where N„ = Max(2N ,i,2jV b + 1) and 5„ = A n / 2 if n is even and S n = 
B( n - 1)/2 if n is odd. However, since we seek a multiplicative splitting of 
(113), a more convenient form to represent Q (A 2 ) is 

s (v) = s„ n (i + (in) 

in which 7 n denote the zeros of the polynomial Si ( — A)*. We immedi- 
ately now identify that each of the product terms in (114) can be factored 
as 

1 + vTEZ d g A/+ (A; 7) M- (A; 7) (115) 

7 
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where 


A+( A; l/ 7 ) = 


s/T^x 


(116) 


M+( A; 7 ) 

is the split function characteristic to the impedance half plane having a 
constant surface impedance I /7 [17]. With the branch choosen so that 
Im ( >/l - A 2 ) < 0, M + (A; 7 ) is explicitly given by 


M+ (A; 7 ) = M- ( — A; 7 ) = < 


A/+ (A; 7 ) Jm( 7)<0 


M+(Ai— y) 


Ln(-y) > 0 , 


(117) 


M+ (cos a; l/?/) 


In this, 


*i(x/ 2 ) [l + y^cos (=^)] [l + v^cos (Sll2=2zi)’ 


y/Z [*, (3tt/2 -a -6)*. (tt/ 2 - a + 6)}‘ 


(118) 


Im(rj) > 0 


cos a 


Im — l/ 7 2 ^ < 0 

# = sin _ 1 (r/) with 0 < Re(0), 


(119) 


and 'I'* (a) is the Maliuzhinets function [18] whose evaluation in algebraic 
form has been given in [19]. 

The determination of Q± (A) is now rather trivial. By substituting (115) 
into (114) we easily obtain 


N s 


Q+ (A) = (?_ (-A) = \fSo M + (A; 7 „) 


(120) 


n= 1 
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Figure 1. (a) Distributed sheet (b) Distributed sheet discontinuity. 



Figure 2. Geometry of the material-material join. 



Figure 3. Illustration of the C contour in the complex X-plane. 




Figure 5. Recessed slab (PEC stub) on a ground plane. 
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Figure 6. H z -polarization backscatter echowidth for a material insert having w=0.04X, 
E=2-j.0001, jj= 1.2 modeled with O(w) low contrast GIBCs. 
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Figure 7. Deformation of the C contour for (a) region 1 integrals and (b) region 2 
integrals. 



Figure 8. Indexing scheme for the constants. 
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Figure 9. H z -polarization backscatter echowidth for material half-plane with t=.20, 
e=2-j.0001, p=1.2. 
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Figure 10. Illustration of a three-layer coated cy linder. 
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Figure 1 1 . Bistatic H-polarization scattering pattern of a circular cylinder of Radius 2.93a. 
coated with a layer 0.07A. thick having e=4 and |i=l; Comparison of fields at a distance 
0.05a. away from the coatings surface. 




Geometry, x > 0 



Tabic 2: Even symmetry parameters associated with the right hand side material x > 0. 
Appendix for the definition of the A/+ split functions. 


































